#include "bernsteinpolyline.h"
#include <qmath.h>
#include "matrix.h"

int fat(int n)
{
    if(n == 0)
        return 1;

    int v = n;
    for(int i = 2; i < n; ++i)
        v *= i;

    return v;
}

int comb(int n, int i)
{
    return (double)(fat(n))/(fat(i)*fat(n-i));
}

BernsteinPolyLine::BernsteinPolyLine(QObject *parent) :
    QObject(parent), m(1,1)
{
}

bool BernsteinPolyLine::bindControlPoints(const QVector<QVector3D>& c)
{
    calcStep(c);
    Matrix contrains(c.size(),c.size());
    int n = c.size() -1;

    for(int i = 0; i <= n; ++i)
    {
        for(int j = 0; j <= i; ++j)
        {
            int l = i-j;
            int c = n-i;

            int sig = j%2==0 ? 1 : -1;

            contrains(l,c) = comb(n,i)*comb(i,j)*sig;
        }
    }


    Matrix points(c.size(),3);
    for(int i =0; i < c.size(); ++i)
    {
        points(i,0) = c[i].x();
        points(i,1) = c[i].y();
        points(i,2) = c[i].z();
    }

    m = contrains*points;

    return true;
}

QVector<QVector3D> BernsteinPolyLine::getPolyLine()
{
    QVector<QVector3D> ret;
    for(float t = 0.0; t < 1.00001; t += step)
    {
        Matrix tt(1,m.nRows());
        tt(0,m.nRows()-1) = 1.0;
        for(int i = m.nRows()-2; i >= 0; --i)
            tt(0,i) = tt(0,i+1)*t;

        tt = tt*m;

        QVector3D v(tt(0,0), tt(0,1), tt(0,2));
        ret.push_back(v);
    }
    return ret;
}

void BernsteinPolyLine::calcStep(const QVector<QVector3D>& c)
{
    float v = 0.0;
    for(int i = 0; i < c.size()-1; ++i)
    {
        v += (c[i]-c[i+1]).length();
    }
    int v2 = qFloor(v);
    v2 = v2 < 1.0 ? 1.0 : v2;
    step = 10.0/v2;
}

QVector<QVector3D> BernsteinPolyLine::getPolyLine(const QVector<QVector3D>& c)
{
    calcStep(c);
    int n = c.size()-1;

    QVector<QVector3D> ret;
    for(float t = 0.0; t < 1.00001; t += step)
    {
        QVector3D v(0,0,0);
        for(int i = 0; i <= n; ++i)
        {
            v += (comb(n, i) * pow(1.0-t, n-i) * pow(t,i)) * c[i];
        }
        ret.append(v);
    }
    return ret;
}
